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Abstract: A fourth-order Runge-Kutta in the interaction pic- 
ture (RK4IP) method is presented for solving the coupled nonlinear 
Schrodinger equation (CNLSE) that governs the light propagation in 
optical fibers with randomly varying birefringence. The computational 
error of RK4IP is caused by the fourth-order Runge-Kutta algorithm, 
better than the split-step approximation limited by the step size. As a 
result, the step size of RK4IP can have the same order of magnitude as 
the dispersion length and/or the nonlinear length of the fiber, provided 
the birefringence effect is small. For communication fibers with ran- 
dom birefringence, the step size of RK4IP can be orders of magnitude 
larger than the correlation length and the beating length of the fibers, 
depending on the interaction between linear and nonlinear effects. Our 
approach can be applied to the fibers having the general form of local 
birefringence and treat the Kerr nonlinearity without approximation. 
For the systems with realistic parameters, the RK4IP results are consis- 
tent with those using Manakov-PMD approximation [11 [21 [3] . However, 
increased interaction between the linear and nonlinear terms in CNLSE 
leads to increased discrepancy between RK4IP and Manakov-PMD ap- 
proximation. 
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1 Introduction 

In optical conimunication fibers, polarization-mode dispersion (PMD, i.e., group bire- 
fringence), chromatic dispersion (CD), and Kerr nonlinearity are the well known eff'ects 
causing signal distortion. In linear region (i. e., without Kerr nonlinearity), the signal 
distortion due to PMD and CD can be predicted and compensated. On the other hand, 
in the limit of pure nonlinear effect (i.e., without linear effects such as PMD and CD), 
the signal amplitude will not be affected and the pulse shape will be unchanged along 
the fiber. Unfortunately both linear and nonlinear effects in real fibers are not negligible. 
To accurately evaluate the signal distortion, one must treat them simultaneously [I]-IS] 
by solving the coupled nonlinear Schrodinger equation (CNLSE), which was proposed 
in 1980's [5] and further studied in many publications. 

For a non-soliton system, the CNLSE is solved with the step size determined by 
dispersion length Ld, nonlinear length L^r, as well as the birefringence related parame- 
ters such as fiber correlation length (Lcorr) (the length within which birefringence axes 
are randomly reoriented) , beating length Abeat (the length determined by birefringence 
strength), and the PMD parameter Dpmd (ps/ Vkm) [l]-[3]. Ignoring the birefringence 
effect, the two dimensional (2D) CNLSE can be reduced to one dimensional nonlinear 
Schrodinger equation (ID NLSE), with its step size being only related with Lu and 
Ln. To efficiently solve the 2D CNLSE and ID NLSE, various approaches have been 
proposed (cf. Refs. [I]-[in| and the references therein). 

In optical communication fibers, the values of A^eat (10^100 m) and L^orr (0.3~^300 
m) are much smaller than Ld and Ljv (up to hundreds of kilometers) PQUJIS]. Dealing 
with the rapidly and randomly varying birefringence in the multispan communication 
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fibers, decreasing the computational step size to a very small percentage of Lcorr and 
Aheat not only needs too much computation but also cannot ensure good enough accu- 
racy because of the computational error accumulated in all steps. To efficiently solve 
the CNLSE, one needs an accurate approach (or algorithm) with large enough step size. 

In Ref. |2] , a coordinate system rotating with the principal axes in each wave plate 
was introduced to get the analytical solutions for the linear and nonlinear effects in 
CNLSE. As a result, the computational step size using the approach of Ref. [2] can be 
increased to a scale significantly larger than Afteat and Lcorr- (Detailed value of its step 
size depends on the interaction between linear and nonlinear effects.) Requirements for 
this approach are that the circular component of the local birefringence in the fiber 
is negligible and that the nonlinear Kerr effect can be approximated as its statistical 
average over the Poincarc sphere (named Manakov-PMD approximation or M-PMD 
approx in this work). Like many approaches used to solve 2D CNLSE or ID NLSE, 
another essential feature of Ref. [5] is that it needs split-step Fourier method (SSFM) , 
which is based on the first order approximation of Baker-Hausdorff formula |16j 

e^e^ = eA+B+k[A,B]+^[A-B,[A,B\\+- ^ (1) 

where [A, B] = AB — BA. Because of this, even in the case of zero birefringence, the 
step size using approach [2] is restricted to 10% ^ 15% of Ld and/or Ljv, assuming the 
computational error tolerance is less than 1% of the pulse peak (cf. the results in l4.2p . 

The concept of interaction picture (IP) was originally used in quantum mechanics. 
Combined with the fourth-order Runge-Kutta (RK4) algorithm, the method of RK4 in 
IP (RK4IP) was recently proposed to study the supercontinuum generation in optical 
fibers (for ID field cases) [inj. Since the IP method does not introduce the error inherent 
in various SSFMs, the computational error of RK4IP is determined by the accuracy of 
RK4. 

The approach of Ref. [TS] cannot be applied directly to the cases with random bire- 
fringence, because in Ref. |15| the linear dispersion operator in ID NLSE was required 
to be constant along the fiber. Motivated by the analytical solution for the linear terms 
in CNLSE, which was proposed in Ref. [2] and is further discussed in the Appendix for 
our calculation, we extend the RK4IP approach from ID NLSE to 2D CNLSE. 

As there is no need to rotate the coordinate system, our approach can be applied to 
the fibers having the general form of local birefringence. Moreover, it treats the Kerr 
nonlinearity without approximation. With the help of the local error method proposed 
in [12], the local error of our RK4IP can be improved to ^ 0{h^), rather than ^ 0{h^) 
of ID RK4IP |15] . As results, for the communication fibers with random birefringence, 
the step size using RK4IP can be orders of magnitude larger than Lcorr and Atieat- 
Without birefringence, the step size of RK4IP can be the same order of magnitude as 
Lo and/or ijv, or, 6 ~ 10 times the step size using the approach of Ref. [2- 

In the parameter regime where communication fibers operate, our results are con- 
sistent with those using M-PMD approx, which was used in many publications (cf. e.g., 
[Hill IS])- Increased interaction between the linear and nonlinear terms in CNLSE will 
lead to increased discrepancy between RK4IP and M-PMD approx. 

2 From CNLSE to M-PMD approx 

2.1 CNLSE expressed in different forms 

To deal with the birefringence related problem, it is convenient to represent a 2D optical 
field with Dirac bra or ket notation |17| . Namely, given a field with its x — y components 
being and Uy, it can be denoted as \u) = {u^,Uy)'^ [or {u\ = {u*,u*)]. Thus, the 
CNLSE discussed in [U [21 13] can be written in the following form with retarded time 
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t = tlab ~ Pt^z: 



j|u;^+2.i A/?|u)+jA/3„|u)t 
or 

j»,+E[A/3|M)+jA/3^|u)t 
where \u) 



/3. 



\u)tt+l 



5 1 1 " 

-{u\u)\u) + -az\u){u\a^\u) + -\v) 



-0, (2) 



{u\u) \u) - -CT31 \u) (u|cri3 |u) 



=0, 



(3) 



dA/3/duj) is related to A, 



beat 



d\u)/dz, \v) = {ulul, ululY, and A/3 (A/3, 

(Lcorr) with A/3 = n/Kbeat {Afiui D p M D / \/^L corr) ■, respectively g]. Here Dpmd is 
the PMD coefBcient (ps/Vkm). The average DGD of a L- km- long fiber can be obtained 
using DGDaug = Dpmd"/L (ps). Obviously, the second term on the left-hand side of 
Eq. ([3]) represents the phase birefringence, while the third term relates to the group 
birefringence (or linear PMD). 

In this work, the three components of the Pauli spin matrices are denoted as 



CTi 



1 

1 



-J 
J 



1 

-1 



(4) 



Generally, the birefringence-induced matrix S in ([2]) and ([3]) has the form of [3] 

E = /3o(z) • (T = ( „ „ ^^"l^'^^ I =(73 cos 61 -Hcti sin 61 cos (/)-f (72 sin 6* sin (/) (5) 
\Pi+JP2 -P3 / 

with (3i = (/3o(2:)|(7j|/3o(2;)) (i = 1,2,3) and \f3Q{z)) [^^{z)] being, respectively, the unit 
vector representing the fiber birefringence in 2D Jones (3D Stokes) space [T7]. For the 
fiber with circular birefringence (e.g., spun optical fiber), /32 7^ in ([S]). 

Parameter fi^^^ = d'^ji/dijj'^ corresponds to the CD parameter at wavelength A by 
/3,,=-CD(A)AV(27rc) (ps/nm-km) with c=3xlOSm/s. 



In Eq. (731 = -(713 = 0-3(71 = ^ ^ 
directly that for a unitary transformation, e.g., 

r 



— j<^2- They are introduced to show 



T 



hi 



12 



tl2 
t*ll 



T- 



'-11 
12 



-tl2 
til 



|tll|' + |tl2p = l, 



(6) 



we have TWisT = (713, provided T ~T* . This means the nonlinear term in the CNLSE 
is covariant for any real rotation. In the following sections, further discussions on CNLSE 
are based on the form of ([3]). 

When the amplitude of \u) is viewed as the square root of optical pulse power, the 
nonlinear coefficient 7 in Eqs. (H])-® relates to Kerr coefficient ri,2, effective mode area 
Aeff, and wavenumber k = 27r/A by 7 = n2k/Aeff. Typical values of A = 1550 nm and 
712 ~ 2.6 X lO^-^^m^/W are used throughout of the paper, unless otherwise noted. 

2.2 M-PMD approx 

Eq. ([3]) can also be expressed as 



Al3\u)+jAp^\u)t 



— mtt+l-{u\u}\u} 



-(m|m)|u)-(73i|u)(u|(7i3|u) 



(7) 



which was named Manakov-PMD equation in Refs. [2J|3]. The last term on the left-hand 
side of ([7]) is the Kerr nonlinearity averaged over the Poincare sphere with the 8/9 factor 
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[U [3J [TH] . The right-hand side of ([7]) is the nonhnear PMD due to ineomplete mixing on 
the Poincarc sphere [U |3] . PhysicaUy it corresponds the rapidly varying fluctuations as 
the polarization state changes [5]. 

In this work, the Manakov-PMD approximation (M-PMD approx) means that the 
variation of the second term on the right-hand side of ^ is approximated by its sta- 
tistical average over the Poincarc sphere [the first term on the right-hand side of ([7])]. 
Thus, the right-hand side of ([7]) is approximated as zero, yielding 



jm 



Al3\u)+jApju), -^\u)u+jI{u\u)\u)=0. (8) 



Without linear birefringence, M-PMD approx ([8]) reduces to Manakov equation p1[2l[T8]. 

The form of Eq. ^ is different from, but equivalent to, the M-PMD approx proposed 

in Refs. [TJ[51|3]. In fact one can obtain the published M-PMD approx, e.g., Eq. (68) of 

Ref. [3], by substituting the unitary transformations \u) = RT\^) into ([5]), with T being 

, , „ / cos a — sin a \ .... „ / — sin a — cos a 

given by mn and R ~ [ . , yieldmg Rz ~ az [ 

° ■' ^ \^ sm a cos a J \ cos a — sm a 

jR^Rz = az02, and [aja. + azAj3]T + jT, = 0. 

3 RK4IP: extension from ID to 2D 

3.1 RK4IP solution for ID NLSE 

Without birefringence, Eq. can be reduced to ID NLSE 

.S-%^g^ + 7H^~ = 0, (9) 

which can be formally viewed as 

3^^+{b + N)u = Q, D = -^^, N = j\u\'. (10) 
Introducing transformation u ~ cxp[j{z — zq)D]u^ , NLSE in the IP has the form [15j 

j— + iVV = (11) 

or 

^=jN'u'^f{zy) (12) 

with 

= exp[-j(^ - zo)^]^cxp[j(z - zo)D] = ij^NIo (13) 

the nonlinear operator represented in the IP. Given u{zn, t), the next step field u{zn+i,t) 
{zn+i = Zn + h) governed by Eq. (fTTj) can be obtained using RK4, with the local error 
of 0{h^). Choosing zq = Zn + h/2 can reduce the required FFTs by 50% (compared to 
the choice of zq = Zn)^ yielding [15] 

u{zn + h, t) ~ exp[j^i)] [m^ -1- fci/i/6 -I- k2h/3 + ksh/S] + k^h/Q 

«n = cxp[j^I)]w(z„,t) 

^'1 = fiZn,ui) = j exp[j^D]N(u{Zn,t))u{Zn,t) 

h = fizn + ^,ui+'^ki) = jN{ul + ^ki)[ul + ^fci] 
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k^ = exp[-j'^D]f{zn+h, ul+hk3)=jN(cxp[j'^D] K + Zi/cg]) cxp[j^L>] K + Zifcs] (14) 



Note that in ([T3|) and (|T4|) . the hncar operator D for ID case was assumed to be 
unchanged in z direction. For a fiber with random birefringence, the transformation 
zq) = exp[j{z— zo)D] = exp[—j{z — zo)l3ujuj{d^/dt'^)/2] in needs to be modified 
correspondly, which is the key point to extend the RK4IP of 15 to 2D case. 



3.2 RK4IP solution for 2D CNLSE with random birefringence 
Eq. ([3]) can also be viewed as 

^2==/3oW-o^(A/?+jm|)-%^^, N2^l[{u\u)+^-aMa,M]- (15) 

Introducing the transformation |m„) = /(z„, zo)|u^), with the operator /(z, zo) given by 
(|?7)) . the CNLSE US]) can be expressed in IP as 



d\u'{z)) 
dz 



jNi{z)\u\z)) = |/(z,u^)), Ni{z)^P{z,zr>)N2{z)i{z,z^). 



(16) 



Given |m„) (the 2D field at z„), the next step field \un+i) can be obtained from Eq. 
by RK4. As in ID case of Ref. [T5], one can choose zq = z„ + /i/2 to minimize the 
required number of FFTs, which leads to 



|u„+i)-/(z„+i,Zo)K+i) = d{^)M2[\ui) + |fcl)^ + |/C2)^ + \k-i)^ 

\ui) = /t(z„,zo)|w„) = dt(_^)M^%„) = d(^)Mi|M„) 
|fci) = |/(z„,<J) 
= j7^^'(^)Mi 

|fc2)-|/(^-„ + ^,uf. + ^fci)) 



(-"„!-(/„) + — ('U„|cri3|u„)J |u„) 



hk-2 



hk2 , 
2 ' 



0-13 , / 



^13/ / 



hki J hki^ 

— kl3|u„ + — , 



hko 



hk-2 



— Fi3|u„+ 2 



1^4 



hki , 
~2~' 



2 ' 



|fc4) = d(2)^^4l/(^« + Kui + hk^)) 



\V^) =d{-)M2\ui + hk3), (17) 



where + cfc^) = + c\ki) [i = 1, 2, 3, c~ -I, /i), 



Ml = mNf,^^rhN^,2-i • • • ^2 • mi, 



M2 = niN^ ■ rriN^-i ■ ■ ■ ri%N^/^+2 ■ 
Ml = mi(-) • m2(-) •• - mAT (-) = rh\ • • • • mj^ (18) 
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and d{h/2) can be obtained according to (P7|) . Introduced in (P7)) . Nh in (|18p is the 
index of the last plate, whereas the unitary matrix rhi in (|18|) and its Hcrmitian (i.e., its 
self-adjoint matrix) mj satisfy m|= ^tt,"^. To order all rhi in one direction, the indexes 
in Ml [given in (jlSp ] does not follow the ordering rule introduced in the Appendix. The 
minus sign in mi(— ) {i = 1,2...) is used to denote 6i < 0, where \6i\ is the length of 
the zth plate discussed in the Appendix. Thus we have rhi{—) = ml. Obviously, when 
7 = 0, (dH) yields |m„+i) = d(|)M2d(|)Mi|u„), which is consistent with ^ [or ^ 
and (j24p ] given in the Appendix. 

4 Applications and discussions 

In the following numerical calculations, the input optical field Win(i) = |win(i)| is as- 
sumed to be a periodic repetition of A'^-bit (N=16) de Bruijn sequence, i.e., Uin(i) = 
J2^=-ac}^B{t — nNT\y), where dB{t) = ^^^q^ aip{t — iTb) and Tb is the time inter- 
val of each bit. Here p{t) determines the elementary input pulse shape and Ui is the 
logic value of the ith bit. Within the time interval [0, Tb], the elementary forms of 
RZ and NRZ pulses (or Marks) are assumed to be p{t) = A/2£'b/2b cos[^ cos^(^)] and 

p{t) = \J Eh/Th^ respectively, with Eh the optical energy per transmitted bit [El [20] . 
Outside this time interval, p{t) is zero. Obviously, Et/Tb is the mark power. To give the 
NRZ optical pulses slightly rounded edges, the input pulses are generated by passing 
through an input optical filter [5]. In this work, it is assumed to be fifth-order Bessel 
type with bandwidth B. 

In all simulations, the fiber loss has been neglected, implying that there are no optical 
amplifiers in the system and, consequently, no amplifier (ASE) noise. 

Only OOK format is considered in this section. Before photodetection, the optical 
signal is assumed to be filtered by a Fabry-Pcrot type channel filter with bandwidth 
Bo = 6.9/Ti, [20l [21]. The square-law-detected signal is then electrically filtered by a 
fifth-order Bessel filter with bandwidth B,. = O.S/Tf, [S] [H] [22] . In the following figures, 
the electrically filtered photoelectric pulses are represented in units of W, since detected 
current corresponds to optical power [2]. 

In our computation, the approach of adaptive step size is used. Namely, given step 
size h, one can use RK4IP (|T7)) to obtain \u{xn+h/2)) / (the fine solution at Xn+h/2) and 
\u{xn + h))c (the coarse solution at Xn + h) from \u{xn))ob (the obtained solution at a;„). 
Similarly, based on \u{xn + h/2)) y, the second fine solution \u{xn + h)) f can be obtained. 
Then the next step size can be adjusted, depending on the difference between |u(x„+/i)) / 
and \u{xn + h))c- According to Ref. |12|. the local error of RK4IP can be improved from 
0{h^) [m to - 0(/i6) by using |m(x„ + h))ob = [16|m(x„ + h)) f - \u{xn + h)),]/l5. 

4.1 RK4IP and SSFM comparison: zero birefringence 

Without birefringence, 2D CNLSE can be reduced to ID NLSE dU) with its split-step 
solution being [h = Zn+i — Zn) 

«(z„+i,t)«F-Me^'^"'''C/], u = u{z,,,t)e^''\<'-^\"\ (19) 

where U is the Fourier transformation of u, while denotes inverse Fourier transfor- 
mation. Note that, in the case of zero birefringence, the SSFM result of Ref. [5] can be 
reduced to 

Obviously, (|19p should yield the same result as the RK4IP solution (ITTI) (without 
birefringence). This is numerically confirmed by considering a NRZ-OOK pulse train 
launched into the fiber with CD=17 ps/(nm-km) and Ae//=80/xm^ [7=1.26 (W-km)"-'^]. 
The bandwidths of the electrical filter and the input optical filter are B^ = 4G and 
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Fig. 1. Electrically filtered pulse as a function of time (T;, = 200ps), with CD=17 
ps/nm-km and fiber birefringence being neglected. The nonlinear coefficient 7 = 



1.26/(W-km) is obtained using effective mode area A^j j 



joa . The input 



mark power E^,|T^, is 20 mW for L=100 km (a) and 2 mW for L=500 km (b). There 
is no visible difference between RK41P solution l|17|l and SSFM solution l|19|l . 



B = 1.75Be. The input mark power is 20 mW for L=100 km (a) and 2 mW for L=500 
km (b). As shown in Fig.[TJ the results of RK4IP without birefringence agrees very well 
with those using SSFM 

4.2 Step size of RK4IP 

Here, the step size means that, within given computational error tolerance (< 1% of the 
pulse peak), the maximum allowable step size at the end of the fiber. 



Table 1. Step size Az using RK41P and SSFM of Ref. [2] obtained with given Ljv and fiber length L. 
The dispersion length 240 km and the computational accuracy < 1% of the pulse peak. 





random birefringence 
RK4IP 


no birefringence 
RK4IP SSFM of |2] 


Az (La, - 40km, L=100 km) 


7 km 


30 km 3 km 


Az (Lat - 400km, L=500 km) 


45 km 


250kni 40 km 



Table 1 shows the step sizes using RK4IP for the fibers of Fig. [T] (a) and (b) with ran- 
dom birefringence [Abeat =50m, Lcorr=10m, and -DpA/£)=1.0ps/(km)^/^, yielding aver- 
age DGD«22ps], compared with the step size using RK4IP (without birefringence) and 
the step size using (|19p . which is the SSFM of Ref. [2] in the case of zero birefringence. 
Notice that the fiber dispersion length Lo ^ l/{(3i^i^Af'^) (A/ is the signal bandwidth) 
for the two cases of Fig. [T]is 240 km, while the nonlinear length ijv ^ l/(7-P) is ~ 
40 km for the case of Fig. [T] (a) and 400 km for (b). 

In each case of zero birefringence, the RK4IP step size is around the smaller one of 
Ld and Ltv, or, about 6~10 times of the step size of SSFM of Ref. [2], which is around 
10% ~ 15% of Lr, and/or Lat. 

Taking into account random birefringence, the RK4IP step size is decreased from 
30km to 7 km for the case of Lat ~ 40 km (L=100km) and from 250 km to 45km for 
Lat ~ 400 km (L=500km). In general, the stronger the interaction between the linear 
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Fig. 2. Tlic filtered photoelectric current (W) vs time (T(, = 100 ps) in a NRZ-OOK 
system consisting of 5 100-km spans of transmission fiber [CDioo=17ps/(nm km), 
A^ff =80 /xm^]. Each span is followed by a 13-km dispersion compensation fiber 
[CDi3=-120ps/(nm-km), A^ff=30 fira^]. Before the 1st 100-km fiber, there is a 
precompensation of -120x6 ps/nm, whereas the 5th 100-km fiber is followed by the 
compensation of -120x5 ps/nm. The input mark power is 10 mW. Birefringence 
parameters are DpMD=2.0 ps/(km)^/^, Lcorr = 10 m, and A^eai = 50 m. There 
is very little difference between the RK4IP solution I I17I I and the M-PMD approx 
using JS}. 



and nonlinear parts in CNLSE, the more impact of Lcorr ^-^id Abeat on the step size, 
assuming that Lo and Ln are much larger than the former. (For the same reason, the 
step size of approach [2] also needs to be decreased correspondingly.) 

4.3 RK4IP solutions and M-PMD approximations with random birefringence 

Fig. [2] shows good agreement between RK4IP result and M-PMD approx for a NRZ- 
OOK system with 5 100-km spans of transmission fiber and 5 13-km spans of dispersion 
compensation fiber. To get the received photoelectric pulses shown in Fig. [21 a pre- 
compensation fiber (6 km) is inserted before the first 100-km fiber, whereas the last 
compensation fiber is 5 km long. The birefringence related parameters are A^eat = 50 
m, Lcorr=^0 m, and Dpmd—'^-O ps/(km)^/^, which means average DGD is around 
2.0v^5(100-H 13) - 2 w47 ps with the birefrin gencc direction being randomly changed 
every 10 m. As in I4.1[ the input optical filter with bandwidth of _B = l.75Be is used 
to generate rounded edges for NRZ signal. Other paratemeters given in the caption of 
Fig. [2] are based on the discussion of Ref. [3]. As plotted in Fig. [21 there is no significant 
difference between the RK4IP solution 1^7^ and the M-PMD approx (|S|). To confirm 
that such agreement is not because of the weak enough nonlinearity, one can increase 
the nonlinear coefficient in ([5]) from |7 — 1.12/(W-km) to I7 = 1.26/(W-km) for the 
transmission fibers and from I7 = 2.99/(W-km) to |7 = 3.36/(W-km) for the compen- 
sation fibers. Our result shows that the discrepancy between the two approximations 
can be more than 20% (not plotted in Fig. [2]). 

Fig. [3] shows the RK4IP solution and M-PMD approx for a RZ-OOK system using 
50 100-km spans of eLEAF fiber [21], with effective mode area being Ae//=72 /im^ and 
CD in the range of 3 ~ 6 ps/(km-nm). To get the received photoelectric pulse shown 
in Fig. [31 each 100-km fiber with CDioo=4.5ps/(km-nm) is followed by -429 ps/nm 
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Fig. 3. Tlie filtered photoelectric current (W) vs time (T^ = 100 ps) in a RZ-OOK 
system using 50 100-km eLEAF fibers [average CDioo=4.5ps/(nm-km), A^ff=72 
fim'^], precompensation of -117 ps/nm before tlie first 100-km fiber, -429 ps/nm com- 
pensation per span, and -234 ps/nm compensation after the 50th 100-km fiber. The 
input mark power is 2 mW and the PMD coefficient is _Dpj\/£3=0.40 ps/(km)^''^. 



0.004 



0.002 



50 lOO^km spans, 2mW, 
DpMD-0-4ps/(l<m)"^ 




RZ-OOK 



n 



I 



7=0, CD=0 
RK4IP 
M-PMD Approx 

1 

1 

'1 

i| 
ll 

'I 



10 



time (T^) 



12 



14 



16 



Fig. 4. Same as Fig. [S] except that CDioo=6.0ps/(nm-km). Because of the change 
of CD, each span is followed by a -585 ps/nm dispersion compensation. The last 
100-km fiber is followed by the compensation of -390 ps/nm. 
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compensation. Also, precompcnsation of -117 ps/nm before the first 100-km fiber and 
last compensation of -234 ps/nm after the 50th 100-km fiber are used. The birefringence 
related parameters are Lcorr = 100 m and Kbeat = 50 m. As shown in Fig. |31 the 
agreement between the RK4IP solution and M-PMD approx is not as good as that of 
Fig. [2] Considering that, due to the linear- nonlinear interaction accumulated in such 50 
100-km fibers, any small change in the related parameters will lead to significant change 
in received pulse shape, the M-PMD approx (dot-dashed) in Fig. [3] can be viewed as a 
reasonably good approximation of the CNLSE solution obtained using RK4IP (dashed). 

Fig. m shows that, when the CD parameter of eLEAF fiber is increased to its worst 
case [6.0ps/(km-nm)], the discrepancy between RK4IP and M-PMD approx becomes 
larger than that of Fig. |31 due to the increased interaction between linear and nonlinear 
terms in the CNLSE. 

5 Summary 

Based on the RK4IP method for ID NLSE [15] as well as the analytical solution for the 
linear terms in CNLSE [2], RK4IP method for 2D CNLSE is presented in this work. 
Without rotating the coordinate system for each step, our approach can be applied 
to a fiber with general form of birefringence. Besides, the Kerr nonlinearity in the 
CNLSE is treated without approximation. Since there is no split-step approximation 
for each step and the local error method of Ref. [12] is used, for normal fibers with 
random birefringence, the step size using RK4IP can be orders of magnitude larger 
than Atieat and Lcom depending on the intensity of the linear-nonlinear interaction. 
Without birefringence effect, the RK4IP step size can be increased to the same order of 
magnitude as Ljj and/or Ljv, or, around 6 ~ 10 times the step size using the approach 
of Ref. [2]. 

In the parameter regime where communication fibers normally operate, our results 
are consistent well with the results using M-PMD approx (|8]) [T] jj] [3]. Increased in- 
teraction between the linear and nonlinear terms in the CNLSE will lead to increased 
discrepancy between RK4IP solution and M-PMD approx. 
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Appendix: Optical field in the fiber without nonlinearity 

For a system with 7 = 0, Eq. ^ can be reduced to 

j^Juiz,t))+Mz)-<?{A^+jA(3^^)\u{z,t))-^^\uiz,t))^0 (20) 

In frequency domain, Eq. (|20p can be simplified as 

j-^jUiz,iv))+[Mz)-a{A(3-A/3^uj) + ^uj']\Uiz,oj))^0, (21) 

where \U{z,u!)) = J \u{z,t))e~^'^'^dt. (Note that in Ref. [2], the Fourier transformation 
was defined to be \U{z,uj)) = J \u{z,t))e^'^*dt.) 
Assuming 

\Uiz,oj)) = A/(z,a;)e^"^'-'(---«)|C/(zo,c^)) = Hz, zo)\U{zo,io)) , (22) 
Eq. (|2T]) is equivalent to the following differential equation for the 2x2 matrix M [2]: 
.dM{z,u}) 



dz 



+ /?o(z)-c?(A/?-A/3^w)Af(z,w) = 0. (23) 
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The correct solution of should be in the form of [2] 



M{Z, Uj)=mNi^ ■ TOJV^-I • • • 7712 • TOl (24) 

m, = eMM^-^-^) (^^i,...,7V,„ 7j = Al3~Ap^uj) (25) 

where 6i = Zi-i — Zi. When Zi > (5^ is the length of the ith plate with its bire- 

fringence direction being represented by Poi^i)- In (pi]) . is the number of the plates 
between zq and z = zq + h and the index Nh denotes the Nhth plate. The expression 
of M{z,uj) given by ((24)) can also be extended to the case of z — zq < 0, provided to 
notice that 1) for each plate, Si < 0; 2) mi (z = 1, • • • ,Nh) are ordered from zq to z, 
i.e., mi is the first left plate of zq. [ In ([25]) rj is determined by the strengths of the 
phase birefringence (A/3) and the group birefringence (A/3tj), which is assumed to be 
independent of z in this work.] 

Using a fixed coordinate system in our approach, the result of (j24p can be simply 
obtained by algebra operations such as [17] 

mi = cos((5i77) + jsin((5i77)/3o(zi) • a = h^q^^ +jfl''^^ ■ cr, 

mim^ = -M« -r^^] [a^^M® - (m» x ^P))] ^ + .7^ • M^^^^ • (26) 

Obviously, in time domain, the transform matrices in (j22p have the form of [i ~ 

/(z, zo) = M(z, f)d(z - zo) 
d(2 - zo)U=2o+'» = er^^^'^ 

M{z,t)^mN,-mN,-i---m2-mi, ^i,^ = e'^^(^^'^^-^) . (27) 
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